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Slowly driven dissipative systems may evolve to a critical state where long pe- 
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(T) , riods of apparent equilibrium are punctuated by intermittent avalanches of activity. 

o 

We present a self-organized critical model of punctuated equilibrium behavior in the 

ON 



context of biological evolution, and solve it in the limit that the number of indepen- 
dent traits for each species diverges. We derive an exact equation of motion for the 



T3 , 

G , avalanche dynamics from the microscopic rules. In the continuum limit, avalanches 

propagate via a diffusion equation with a nonlocal, history-dependent potential rep- 
resenting memory. This nonlocal potential gives rise to a non-Gaussian (fat) tail for 
the subdiffusive spreading of activity. The probability for the activity to spread be- 
yond a distance r in time s decays as ^/^'s _3//2 x 1 / 3 exp [— la: 1 / 3 ] for x = ^ 3> 1. The 
potential represents a hierarchy of time scales that is dynamically generated by the 
ultrametric structure of avalanches, which can be quantified in terms of "backward" 
avalanches. In addition, a number of other correlation functions characterizing the 
punctuated equilibrium dynamics are determined exactly. 
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I. INTRODUCTION 



Many natural phenomena evolve intermittently rather than following a uniform, gradual 
path. In particular the dynamics out of equilibrium may follow a step-like pattern with long, 
dormant plateaus interrupted by sudden bursts, or avalanches, where the accumulated stress 
is released. Avalanche dynamics violates the picture of gradualism where large systems 
evolve continuously, for instance to a local energy minimum. The bursts which separate 
subsequent metastable states may eventually span all scales up to the system size |T| thus 
providing a general mechanism for long range spatiotemporal correlations. The appearance 
of fractals, 1/f noise, Levy flights, etc., have been unified by relating these phenomena, in 
a broad class of non-equilibrium models, to an underlying avalanche structure f2|. 

Even though the theory of uniformitarianism, or gradualism, has historically dominated 
both geology and paleontology, prototypical examples of avalanche dynamics lie in these 
two domains 0. For instance, the distribution of earthquake magnitudes follows a power 
law found by Gutenberg and Richter [|]]. The scale- free variation from small events to 
large events indicates a common dynamical origin and eventually led to the suggestion that 
earthquakes are an example of self-organized criticality ||. Most of the time, the crust of the 
earth appears stable. These periods of apparent equilibrium are punctuated by earthquakes, 
which take place on a fractal fault structure |]] that stores information about the history of 
the system. 

Actually, over twenty years ago, Gould and Eldredge proposed that biological evolution 
also takes place in terms of punctuations, where many species become extinct and new 
species emerge, interrupting periods of low activity J7|. In this context, punctuated equi- 
librium usually refers to the intermittent dynamics of single species, where morphological 
change is concentrated in short intervals in time interrupting long periods of stasis. These 
punctuations may be correlated to large extinction events in the global ecology, which may 
themselves be distributed according to a power law analogous to the Gutenburg-Richter law 
for earthquakes 



2 



This view was promoted by Bak and Sneppen [[| who introduced a simple self-organized 
critical (SOC) model for coevolutionary avalanches of different species in an ecology. The 
model explicitly treats macroevolution as a many-body statistical problem where the fitness 
landscape in which each species evolves is affected by changes in other species in the ecology 
The mutation of the "least fit" species in the system completely changes the fitness landscape 
of some other species and coevolution with punctuated equilibrium behavior is obtained. 
This occurs without fine tuning parameters and without the need for external shocks leading 
to cataclysms of mass extinction. The extinction events (avalanches) have a power law 
distribution of sizes, where most extinction is concentrated in the largest events. This 
provides some theoretical underpinning for catastrophism rather than uniformitarianism 
as the defining characterization of evolutionary history ||10|| . The model may capture some 
robust features of real evolution, such as punctuated equilibrium and catastrophism, in spite 
of its drastic simplifications. 

In fact Ito |Tl|] has related the spatiotemporal pattern of earthquakes in California to the 
avalanche pattern found in the Bak-Sneppen evolution model. The data for the pattern of 
successive earthquakes over time may be fractal in space and time, so that each earthquake 
can be viewed a single event within a much larger avalanche structure consisting of many 
earthquakes. If this picture is correct, earthquakes can be described as a "fractal renewal 
process" p|,|l"2f with a power law distribution of times between subsequent earthquakes in a 
given region. The analogy can be made by replacing the least fit species with the weakest 
link in the earth's crust. Ito's observations of return times for earthquakes in California are 
consistent with general scaling relations valid for the Bak-Sneppen model, as well as other 
extremal SOC models 0] for invasion percolation fl3 |, interface depinning [14]], and flux 
creep [fT5] . 

Although the Bak-Sneppen model is an extremely simple model of SOC and a few exact 
results as well as many scaling relations are known for it [FJ], it has not (yet) been solved. A 
variety of similar evolution models have been subsequently introduced to incorporate more 
aspects of reality such as speciation and external shocks |16 . These models have been 



primarily studied numerically. 

Recently, we have proposed a multi-trait evolution model with M independent, internal 



degrees of freedom for each species []T7|]. This model is solvable in the limit M — > oo; it 
includes the Bak-Sneppen model when M — 1. We have previously reported some of the 
simplest results [fL7| . Here we present a more complete analysis of the model. Note that 
the M — > oo limit is not mean-field theory |18]| because it contains spatial correlations and 



punctuated equilibrium dynamics. In addition, we find rigorous results that might be too 
delicate to detect with numerical calculations alone, such as the non-Gaussian tail for the 
distribution of activity. 

In the multi-trait evolution model, we have incorporated the notion that the survivability 
of each species depends on a number of independent traits associated with the different tasks 
that it has to perform in order to survive [|19| . Evolution proceeds via an extremal dynamics 



where the species in the global ecology with the lowest barrier to mutation, or the least 
fit, mutates. This event affects certain barriers to mutation or fitnesses of other species 
in the system which are related through, for instance, a food chain. As a consequence of 
the interaction between species, even species that possess well-adapted abilities, with high 
barriers, can be undermined in their existence by weak species with which they interact. This 
may lead to a chain reaction of coevolution. The pattern of change of individual species 
exhibits punctuated equilibrium behavior which comes from episodes of mass extinction 
events sweeping through the species. Punctuated equilibrium is described by a Devil's 
staircase, as shown in Fig. [I]. The introduction of many internal traits for each species 
is consistent with paleontological observations indicating that evolution within a species is 
"directed" ||19|| ; morphological change over time is concentrated in a few traits, while most 



others traits of the species are static. 

Our main analytical results are as follows: For M — > oo, we derive an exact equation 
of motion, given in Eq. ( |3.7| ), for the macroscopic observables in the SOC state from the 
microscopic dynamics. From this equation, which is our central result, one can extract 
separate equations for the temporal and spatial distribution of avalanche sizes, both of which 



are determined exactly in Eqs. ( [2.2| , |2.3| ) and were previously derived by us in Ref. [17 



In the continuum limit, the subdiffusive dynamics is given by a "Schrodinger" equation 
with a nonlocal potential in time, Eq. ( |3.8| ). This potential represents memory. The exact 
asymptotic solution, Eq. Q3.9|), of the Schrodinger equation at large length and time scales 
has a non-Gaussian tail. This solution describes the nontrivial spatiotemporal pattern of 
activity in the self-organized critical state. The anomalous diffusion arises from a long-term 
memory effect due to the ultrametric tree structure of avalanches. This tree structure of 
activity is quantified by calculating the ultrametric distances between subsequent extremal 
(minimal) sites. The probability distribution for this distance is a power law at large times 
and asymptotically approaches the probability distribution of all backward avalanches. This 
latter quantity, which we will define in the text, is calculated exactly in Eq. ( |5.3|) . A number 
of other distribution functions are also determined. The critical exponents are D = 4, 
r = 3/2, t r = 3, 7 = 1, v = a = 1/2, rf = 2, rf = 3/2, and r FIRST = 2 - d/4. 
The nonlocal time dependence of the dynamical equations and the ultrametric structure of 
avalanches suggest a possible relation between glassy dynamics and self-organized criticality 

In Section II, we introduce our model and show that it self-organizes to the critical state. 
It is demonstrated that avalanches in the critical state have an ultrametric tree structure. 
In Section III, we derive the equation of motion for the critical state for the M — > oo model 
and present our main analytical results for the anomalous diffusion. In Section IV, we show 
that the critical behavior is characterized by simple power laws with specific exponents that 
verify general scaling relations for nonequilibrium phenomena. In Section V we calculate 
the distributions for "backward avalanches." Due to irreversibility these differ from the 
usual "forward" avalanche distributions. The backward avalanches can be related to the 
ultrametric distances between subsequent activity. Most details of our calculations have 
been deferred to the Appendices. 
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II. THE MULTI-TRAIT EVOLUTION MODEL 



Our model is defined as follows: A species is represented by a single site on a d- 
dimensional lattice. The collection of traits for each species is represented by a set of 
M numbers in the unit interval. A larger number represents a better ability to perform that 
particular task, while smaller numbers pose less of a barrier against mutation. Therefore, 
we "mutate" at every time step the smallest number in the entire system by replacing it by 
a new (possibly smaller) number that is randomly drawn from a flat distribution in the unit 
interval, V. Choosing the smallest random number mimics the Darwinian principle that the 



least fit species mutates p2 



The dynamical impact of this event on neighboring species is simulated by also replacing 
one of the M numbers on each neighboring site with a new random number drawn from V. 
Which one of the M numbers is selected for such an update is determined at random, since 
we assume that a mutation in the traits of one species can lead to an adaptive change in 
any one of the traits of a dependent species. The interaction between the fitnesses of species 
leads to a chain reaction of coevolution. 



A. Self-organization 

The sequence of selective random updates at extremal (minimal) sites with nearest- 
neighbor interactions drives the system from any initial state to a self-organized critical state 
in which species exist in a state of punctuated equilibrium with bursts of evolutionary activity 
that are correlated over all spatial and temporal extents. In this state almost all species have 
reached fitnesses above a SOC threshold, enjoying long periods of quiescence, interrupted 
by intermittent activity when changes in neighboring species force a readjustment in their 
own barriers. A snapshot of the stationary state in a finite one dimensional system is shown 
in Fig. |[ 

The self-organization process for the finite M model is similar to that for the Bak- 
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Sneppen model (M = 1). It is described by a "gap" equation that relates the rate of 
approach towards the stationary attractor to the average avalanche size. This equation 
demonstrates that the stationary state of the system is a critical state for the avalanches, 
where the average avalanche size diverges. Following Ref. |2| we define the gap, G(s), to be 
the largest of the minimal random numbers selected up to time s given an initial state at 
s = where all the random numbers are uniformly distributed in the unit interval. Therefore, 
G(s = 0) = 0{1/ ML d ). As evolution proceeds, the gap G(s) increases monotonically in a 
stepwise fashion with intermediate plateaus that become longer and longer. These plateaus 
occur when subsequent minimal random numbers in the system are smaller than a minimum 
at some previous time. One can assign to each plateau an avalanche which starts and ends 
with consecutive changes in G(s), and which consists of all the random numbers in the gap 
below G(s). As the gap increases, the probability for the new random numbers to fall below 
the current gap increases also, and longer and longer avalanches typically occur. 

Following Ref. |2j again, it can be shown from the law of large numbers that in the limit 
of large system sizes L, the growth of the gap versus time s obeys the gap equation: 

dG(s) 1 - G(s) 



ds ML d < S > G(s) 



(2.1) 



As the gap increases, so does the average avalanche size < S >g{s), which eventually diverges 
as G(s) — > G c . In the limit L — > oo, the density of sites with random numbers less than G c 
vanishes, and the distribution of random numbers is uniform above G c . The gap equation 
( |2.1D contains the essential physics of SOC phenomena. When the average avalanche size 
diverges, < S >g(s)— * oo, the system becomes critical. At the same time, approaches 
zero, which means that the system becomes stationary. For any finite M the process of 
self-organization is the same as for the M — 1 model, and all of the results derived for that 
case apply. Since M just enters as a rescaling of the system size L in Eq. ( |2.1| ), it is plausible 
that the critical behavior for any finite M is in the same universality class as for M — 1. 
We have shown |T7[ that the limit M — > oo is a different universality class. 
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B. Ultrametricity of Avalanches 



It is useful to consider the case where, at a certain time, the smallest random number in 
the system has the value A. A A avalanche by definition consists of all subsequent random 
numbers which are below A. The A avalanche that started at s ends at the first instant 
s + s' when the smallest random number in the system is larger than A. All of the random 
numbers that are below the threshold value A at any one instant in time are called "active" 
because they make up the A avalanche. 

We now consider the sequence of minimal values A m j„(s) comprising any A avalanche. 
Each value A m j n (s) has a parent barrier value A m f„(s — s') preceding it within the A avalanche. 
This parent is the barrier that introduced the particular random number into the system 
that became the minimum at time s. Obviously, the parent of the A m j n (s) value also has 
its own parent. One can therefore place all of the barrier values within a given avalanche 
onto a tree, as shown in Fig. |[ The distance on the tree between any two active barrier 
values at a given time is determined by the most recent common ancestor of the two values. 
This distance is ultrametric f20| . In Section V we relate the probability distribution of the 
ultrametric distances between subsequent minimal random numbers to the distribution of 
all backward avalanches, which is determined exactly. 

For finite M, active barriers can be eliminated both by becoming the global minimum at 
some time, or by being chosen in the nearest-neighbor interaction when the global minimum 
occurs on a neighboring site. The case M — > oo of the model is special because the existing 
active barriers that any species possesses can only be changed if these barriers themselves 
become the global minimum, despite the nearest-neighbor interaction in the model. Since 
there are infinitely many barriers on each site, no existing active barrier is ever likely to be 
chosen for an update in a nearest-neighbor interaction. The nearest-neighbor interaction 
can only create new active barriers. This fact allows us to formulate a cascade mechanism 
which can be solved exactly. 

In what follows we shall consider the stationary state behavior in the limit of infinite 



S 



system size, and confine our presentation to the d = 1 dimensional case of our model 
with M — ► oo. We point out later which critical exponents are dimension dependent and 
dimension independent, and we will discuss in more detail the results for higher dimensional 
lattices elsewhere ||23|| . To simplify the algebra further, we make a slight modification of 
the model without restricting the generality of the results: At each time step during the 
avalanche, the smallest active barrier is set to unity instead of being replaced by a new 
random number. Then, there is either no new active barrier created with probability (1 — A) 2 , 
or one such barrier is created to the left or to the right of the minimum with probability 
A(l — A) in either case, or two active barriers are created to the left and the right of the 
minimum with probability A 2 . 

The spatial and temporal correlations in our model can be separated into two independent 
equations of motion for the width and duration of avalanches. For instance, the distribution 
of avalanches sizes s is given by (see Appendix A) 

P x (s) = 4 s \ s -\l- \) s+1 A(s) 

A(s) = —f-i V— . (2.2) 

r(i)i> + 2) 

For A = A c = 1/2, P\ c (s) ~ s~ T for large sizes with r = 3/2. For A < A c , the critical 
avalanches are subdivided into smaller avalanches and the distribution acquires a cutoff. 
The critical behavior is quite different than the M = 1 model where the numerical result 
is r ~ 1.1 in one dimension [0. The model for M = oo clearly represents a different 
universality class than the original Bak-Sneppen model. Actually, the temporal behavior in 
this case is identical to the temporal behavior of the random-neighbor evolution model [|I8| 



Unlike the random-neighbor model, though, the M — > oo model also exhibits spatial 
correlations, leading to punctuated equilibrium behavior. In Ref. [DJ we found that the 



probability f r of a A c avalanche to ever affect a particular site of distance r from its origin 
is exactly given by 



implying that the probability for an avalanche to reach precisely to a site of distance r falls 
asymptotically as r~ TR with = 3. Noting that one particular site can not fully contain 
the spread of an avalanche, we have obtained an equation for the probability f(0,r) of a 
A c avalanche to start at the origin i = and to ever leave a box of radius i = r around its 
origin. The leading asymptotic behavior of the solution for confined avalanches is given by 



confirming that tr = 3. The calculation leading to Eq. ( p.4|) is given in Appendix B. A 
comparison of numerical results for f(0,r) with the asymptotic behavior, given in Fig. |], 
show perfect agreement. 



Ultimately, spatial and temporal correlations are interrelated through the microscopic 
rules for the propagation of activity in space and time. Ideally, one would like to determine 
the propagator G(r, s) which is the probability that the minimum barrier value will be 
located at position r at time s given that it was at location at time for the infinite A c 
avalanche. We have not been able to analytically calculate this quantity directly. We have 
instead focussed on a different quantity. Let F\(r, s) be the probability for a A avalanche to 
survive precisely s steps and to have affected a particular site of distance r from its origin. 
Conceptually, the quantity F\ c (r, s) may roughly correspond to an envelope function of the 
propagator G(r, s). This relation between F and G is explained in Fig. [5]. Due to the lack of 
any scale in the model, it is plausible that the asymptotic behavior of G and F is identical, 
as comparison with numerical calculations suggests. 

The direct analysis of this envelope function proves to be rewarding in many respects. 
We find a cascade equation which can be reduced to separate equations for spatial and 
temporal correlations. In the continuum limit, the avalanche dynamics is given in terms 
of a Schrodinger equation with a nonlocal potential. We solve this equation to find the 




(2.4) 



III. THE CASCADE EQUATION AND MAIN RESULTS 
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leading asymptotic behavior of F\ c (r,s). This calculation yields a non-Gaussian tail in the 
distribution with an avalanche dimension D = 4, signaling subdiffusive behavior for the 
spread of activity. 

First, we consider P\(r,s), the probability that the A avalanche dies precisely after s 
updates and does not affect a particular site r away from the origin of the avalanche. The 
quantities P\{r, s) and F\(r, s) are simply related: 

P A (r, s) = P A (r = oo, s) - F x (r, s), (3.1) 



where P A (r — oo, s) — P A (s), as given in Eq. (|2.2| ). Since we consider the avalanche to 
start with a single active barrier at r = and s — 0, it is P\(r — 0, s) = for all s > 0, 
and P\(r,s = 0) = for all r. The remaining properties of a A avalanche can be deduced 
from the properties of avalanches that ensue after the first update. It will terminate with 
probability (1 — A) 2 after the first update when the update does not produce any new active 
barriers. Thus, it is P\(r,s = 1) = (1 — A) 2 . For avalanches surviving until s > 2 we find 
for r > 1 

P A (r, s) = A(l - A) [P x (r - 1, s - 1) + P A (r + 1, s - 1)] 

+A 2 Px(r - 1, s')P x (r + 1, s - 1 - s') (3.2) 

s'=0 

in the following way: The first update may create exactly one new active barrier with 
probability A(l — A) either to the left or to the right of the origin (i. e. one step towards 
or away from the chosen site of distance r). In this case, the properties of the original 
avalanche of duration s are related to the properties of an avalanche of duration s — 1 with 
regard to a site of distance r — 1 or r + 1, respectively. Finally, the first update may create 
two new active barriers with probability A 2 to the left and the right of the origin. Then, 
the properties of the original avalanche of duration s are related to the properties of all 
combinations of two avalanches of combined duration s — 1. Both of these avalanches evolve 
in a statistically independent manner for M = oo. Since only one of these avalanches can be 
updated at each time step, their combined duration has to add up to s — 1 for this process 
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to contribute to the avalanche of duration s. For any such combination, the probability to 
not affect the chosen site of distance r from the origin is given simply by the product of the 
probabilities for the two ensuing avalanches to not affect a chosen site of distance r — 1 or 
r + 1, respectively 

Before proceeding with the solution of Eq. fl3.2| ), we review some limiting cases to obtain 
previously derived results jl7| . Considering the limit r = oo, for s > 2, Eq. (|3.2|) reduces to 

P x (s) = 2A(1 - X)P x (s - 1) + A 2 ]T P x (s')P x (s - 1 - s'). (3.3) 

s'=0 

This is the cascade equation for the lifetime distribution of avalanches whose solution is Eq. 
( f2~2|) (see Appendix A). Near the critical point, A A = A c — A — ► + , the lifetime distribution 
obeys a scaling form 

P x (s)~ s ^ G (sA\^ (r = ^,a= 1 -). (3.4) 



Similarly, we can rederive Eq. ( [2.3|) for the spatial correlations. Defining N x (r) 
E s P\{r, s), Eq. (|J) y ields N\(0) = 0, and for r > 1 



N x (r) = (1 - A) 2 + A(l - A) [N x (r - 1) + N x (r + 1)] + \ 2 N x (r - l)N x (r + 1). (3.5) 

N x (r) is the probability for an avalanche of any temporal extent not to reach a particular 
point of distance r before it dies. Eq. ( |3.5|) for A = A c = 1/2 is solved by N Xc (r) = 1 — f r 
with f r given in Eq. (|2.3|) . Close to the critical point this quantity also obeys a scaling form: 



1 - N x (r) ~ -±^H (r A\») (r R = 3,u = ^j. (3.6) 



The equation governing the envelope function F x (r, s) is obtained by inserting Eq. ( |3.1|) 
into Eq. It is F x (0,s) = P x {s), F x (r,0) = 0, F x (r > 1, s = 1) = 0, and for s > 1, 

r > 1, 

F x (r, s + 1) = A(l - A) [F A (r - 1, s) + F A (r + 1, s)\ 

+A 2 ^ P A ( S - s') [F x (r - 1, s') + F x {r + 1, s')\ 



s'=0 
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-A 2 J2 F x(r ~ 1, s')F x (r + 1,8-8'). (3.7) 



s'=0 



Now we will focus on the spatiotemporal correlations at the critical point A = A c = 1/2. 
For sufficiently large values of r and s we will show that F\(r, s) — > for r — > oo sufficiently 
fast such that we can neglect the nonlinear term in Eq. (|3.7|) . We can take the continuum 
limit and obtain |24| 

~ \ V 2 r F(r, s) + \[ V(s - s')F(r, s')ds', (3.8) 
a "Schrodinger" equation in imaginary time, for F(r, s) with a nonlocal memory kernel 



V(s) = P(s) — 25(s), where P(s) = P\ c (s) given in Eq. ( |2.2j ), and 5(s) is the usual Dirac 
8 function. Note that it is the statistical independence of the avalanches that gives V(s) 
in terms of the probability distribution of avalanche sizes. The memory in the system is 
characterized solely in terms of this distribution. 

This nonlocal potential with the integral kernel V(s) contains all of the history depen- 
dence of the process. In its absence the system would have no memory and be purely 
diffusive with a Gaussian tail F ~ e~ r ' 2//2s . In its presence the probability to have reached 
a site at distance r at time s gets contributions from avalanches that reached r at earlier 
times s' < s. These contributions are weighted according to P(s — s') which has a power 
law tail. Avalanches contributing to F(r, s) consist of sub-avalanches, one of which reaches 
r in time s' while the other's combined duration is s — s'. The sub-avalanche tree structure 
gives a hierarchy of time scales. This changes the relaxation dynamics to be non-Gaussian, 
and we find as our main result that 

^'e" 1 ^ (r 4 >s>l) . (3.9) 

It is immediately clear from the form of the solution that the avalanche dimension for the 
subdiffusive spreading of activity {r D ~ s) is D = 4. Diffusion is slowed down because the 
activity has a tendency to revisit sites, and the system remembers these previously visited 
sites. Considering x = r D / s as the scaling variable, the variation of the exponent is much 
slower with a "fat" tail (~ x 1//3 ) compared to a Gaussian tail. 
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In Appendix C we derive the complete leading asymptotic behavior given in Eq. fl3.9|) . 
Here we will just show how the history dependence in the Schrodinger equation ( |3.8| ) gives 
rise to the non-Gaussian tail in the exponential of Eq. ( |3.9| ) . 



A. Calculation of the Fat Tail 



Using a Laplace transform, F(r, y) = J °° ds e ys F(r, s), Eq. (|3.8|) turns into an ordinary 
second-order differential equation in r, 

V r 2 F(r, y) ~ [2y - V(y)] F(r, y), (3.10) 

where V(y) is the Laplace transform of V(s). Since F(r, s) is falling for large r, it is 



F(r, y) ~ C{y) exp 



-ryJ2y-V(y) 



(3.11) 



and we assume that C(y) is a sufficiently well-behaved function near y = 0. 

The inverse Laplace transform yields a contour integral with a contour extending just to 
the right (77 > 0) of the imaginary y-axis: 

>+v dy 



F(r, s) ~ 



d+t? 2m 



C(y) exp 



ys -r<J2y- V(y) 



The limit of large s corresponds to small values of y where V(y) 



(3.12) 



-2^/y such that 



2y — V(y) ~ \[2y x l^ . Note that the contribution from the left-hand side of Eq. ( p.8|) 
does not effect the leading order which merely consists of a balance between the terms on 
the right-hand side. After rescaling y — > y/s it is 

(7(0) 



F (r, s) ~ 



dy exp 



• ( — J // ' 

S4 



(3.13) 



2-Kis Jc 

where C is a small piece of contour that crosses the real y-axis from below just to the right of 
y = 0. It emerges that F(r, s) is exponentially cut off when r 3> s 1//4 3> 1, determining that 
the avalanche dimension is D = 4. In that limit, we can perform a steepest-descent analysis 



of the integral p5| . First, we note that the expression in the exponent has a moving saddle 
point at y = It- 4 / 3 ^ 1 / 3 . To fix the saddle point, we substitute y = yov to obtain 
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F(r,s) ~ I dvexp y (v - 4vt) , (3.14) 



2nis Jc 

where the contour C is deformed such that it crosses the real axis on the saddle point at 
vq = 1. To find the steepest-descent path for the contour in the neighborhood of the saddle 
point, we set w ~ 1 + e + iS with e, 5 <C 1. Thus, in the exponent, it is v — 4-y 1 / 4 ~ 
(—3 + |e 2 — |5 2 ) — % (f e^) , indicating that the steepest-descent path (which is always also 
a constant-phase path) is given by e = in the neighborhood of the saddle point. Thus, we 
substitute v = 1 + iS in that neighborhood and obtain the non-Gaussian tail in Eq. ( ^9|) , 



F(r, s) ~ C'(r, s) exp 



3 /r 4 V 



r 4 > s > 1), (3.15) 



where C'(r, s) only contains powers of r and s. The function C'(r, s), as well as the behavior 
of F(r, s) for s ^> r 4 ^> 1, is determined in Appendix C. In Appendix D we will also consider 
the distribution for avalanches that are fully confined in a box of size r (see Fig. |5|). We find 
that the dominant relaxation behavior is also given by Eq. ( |3.15| ). 



B. Numerical comparison with the propagator 

We are now in a position to attempt to compare these results with numerical measure- 



ments of the actual propagator G(r,s). Based on earlier numerical investigations p6 
assume that G(r, s) has the scaling form 



G(r, s) ~ r d g \- (r*/s > lj . (3.16) 

At each update s we determine the location r of the minimal random number in a surviving 
avalanche that started at s = 0, r = 0. We bin counts as a function of r and r 4 /s for 
avalanches of duration s < 10 8 . In Fig. ^|we plot the numerical results for G(r, s) for various 
values of r A /s as a function of s. This makes the numerically imposed cut-off at s — 10 8 
apparent. The data indicates that G(r, s) rises linearly with r for each value of r A /s; thus 
5=1. In Fig. |7| we plot G(r, s)/r vs. s for different values of r 4 /s for 10 4 < s < 5 x 10 7 . 
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Note that all data "collapses" onto plateaus whose mean value and standard deviation is 
given by the crosses with error bars to the right of each curve. 

We attempted to numerically determine the asymptotic tail of the function g for large 
values of x = r 4 / s by assuming in accordance with Eq. ( |3.9| ) that 

g(x)~ e - Ax \ (3.17) 

From the sequence of plateau values in Fig. [7| for g{x) we determine a sequence of extrapolants 

In (— hiq(x)) In A . 

\ yK >> ~ a+ - , 3.18 

In a; Ini 

In Fig. H we plot this sequence of extrapolants as a function of 1/ In x and estimate that a = 
0.35 ± 0.03, in reasonable agreement with with the value a = 1/3 from the exponential fall- 
off for F(r, s). Thus, the behavior of G(r, s) is consistent with the non-Gaussian asymptotic 
behavior of the envelope function F(r, s). 



IV. SCALING RELATIONS 

In the SOC state, spatial and temporal correlations are interrelated. This interrelation 
is expressed via scaling relations. In a broad class of SOC models, including the evolution 
models, the knowledge of just two scaling coefficients, such as r and D, is sufficient to 
determine any other known coefficient of the SOC state, including the approach to the 
attractor, through these scaling relations ||. For example, we have shown in Sec. Ill that 
the activity in the SOC state spreads in a subdiffusive manner, r ~ where D = 4 is 
the avalanche dimension. In Ref. |17| we have numerically determined the root-mean-square 



distance of the location of the activity at time s and shown that it scales as s 1//4 . In Fig. |9| 
we show that the number of sites covered by the activity also grows as s 1 / 4 . This verifies 
that rather than being multifractal there is only one dimension for the avalanche. 

The probability distribution of A c avalanche duration is a power law with exponent 
r = 3/2, and we have shown that the probability distribution of spatial extents of avalanches 
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is also a power law with critical exponent tr = 3. This verifies the scaling relation tr — 1 = 
D(r — 1). It is easy to show that the average size of an avalanche diverges with exponent 
< s >~ (AA)~ 7 with 7 = 1, and previously we derived that the critical exponent a describing 
the cutoff of avalanche sizes below the critical point is 1/2. These results verify the scaling 
relation 7 = (2 — r)/a. The result in Sec. Ill for the cutoff in the correlation length v = 1/2 
in Eq. ( |3.6|) verifies the scaling relation v = l/(aD). Our analysis of the model in higher 
than one dimension shows that all of these above mentioned exponents are independent of 
dimension. 

The punctuated equilibrium behavior, however, does depend on dimension. Each site 
is visited many times in a long lived avalanche. The intervals between subsequent returns 
to a given site correspond to periods of stasis for a given species. As shown in Fig. [1], the 
accumulated number of returns to a given site forms a "Devil's staircase"; the plateaus in 
the staircase are the periods of stasis for that species. The punctuations, i. e. the times 
when the number of returns increases, occupy a vanishingly small fraction of the time scale 
on which the evolutionary activity proceeds. The distribution of plateau sizes is the same 
as the distribution of first returns of the activity to a given site, -Pfirst(s)- We found in 
dimensions d < 4 that P F irst(s) ~ s - t first f or i arge s w j t h r FIRST = 2 - d/D 0. For 
M — > 00 in d = 1 dimension, the scaling relations therefore predict tfirst = 7/4. We have 
measured tfirst = 1.73 ± 0.05. 

A. Forward Avalanches 

The quantity P\{s) is the conditional probability to have a forward avalanche of size s 
given that the signal at the starting point was equal to A. Such avalanches are defined by 
looking at the first moment forward in time when the signal, \ m in exceeds its current value 
A. It is important to note that this conditional forward probability distribution is exactly 
equal to the punctuating avalanche distribution. Punctuating A avalanches are defined as 
the intervals separating subsequent punctuations of the barrier A by the signal A mi „(s); see 
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Fig. [IJJ. The probability distribution for the A avalanches does not depend on the precise 
value of the site that started the A avalanche (as long as it is > A) because this value is erased 
from the system after the first update. Thus, the minimal numbers selected at different time 



steps are distributed as 



p{Kun > A) = 1/ < S > X 

oo 

<s> x =J2sPx(s) . (4.1) 



s=0 



Substituting Eq. (gj) gives < s > A = 1/(1 - 2A), and p(X m in = A) = 2 for < A < 1/2. 

The distribution of all forward avalanches, Pf l {s) is obtained by integrating the condi- 
tional probability from Eq. (|2.2|) with the proper weight p{\ m in = A) for the starting value 



of the avalanche [27L2S3]. This distribution measures avalanches which begin at every time 
step, s', and end at the first moment forward in time, s' + s, when A miri (s' + s) > A m j n (s'). 
We find 

Pf(s) = ti ,2 p{\ mm = \)P x {s)d\ + ^ , (4.2) 

which agrees with numerical simulations. For large s this distribution is a power law 

Pf{s) ~ s~ T f l where rf = 2 . (4.3) 

This particular exponent rf 1 turns out to be superuniversal and equals 2 for a wide variety 
of extremal models |28]j2||. 



V. BACKWARD AVALANCHES AND ULTRAMETRICITY 

The properties of the system under time reversal can be studied in terms of 
backward avalanches. Now we look for the first moment back in time when 
X m i n (s' — s) > \ m i n (s') = A ]27LpS],B]. The definitions of forward and backward avalanches 
are illustrated in Figs. |TD| and [IT]. These figures demonstrate the hierarchy in the avalanche 
structure: all forward and backward avalanches that start inside one big forward (backward) 
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avalanche are constrained to not go beyond the limits of the parental avalanche and, there- 
fore, can be considered to be its sub-avalanches. Each sub-avalanche in turn has its own 
sub- avalanches, and so on. One can look at the entire activity as one great parental critical 
avalanche, which began in the distant past. It contains sub-avalanches of all sizes. 

We can determine the distributions for A backward avalanches -P^(s) and all backward 
avalanches P£ ll (s) exactly. Suppose we have a temporal sequence X m i n (s') which is an 
ensemble of N, A punctuating avalanches, where N is a sufficiently large number. The 
average number of A punctuating avalanches of size s in such an ensemble is given by 
N(s) = NP\(s). At the end of any A avalanche of size s, precisely (s + 1) new random 
numbers have been introduced into the system that were not there when the avalanche 
began. All these random numbers are uncorrelated and uniformly distributed between A 



To have X m i n (s' + s) = X we need the minimal number in the system to lie between A and 
A + dX. This number can be only in this set of new random numbers, since at the beginning 
of the A avalanche every number in the system was by definition larger than or equal to 
A. The probability that at least one of these new numbers will be between A and A + dX 



will be changed by dN = -N(s + 1)^. Of course the sum of the temporal durations of 
these avalanches will remain constant. This leads to the following differential equation for 
the average size of an avalanche 



which is analogous to Eq. (16) in Ref ||, and also shows that the avalanche size diverges 
with exponent 7 = 1. 

The number N^(s) of valid A backward avalanches of size s in our ensemble of punc- 
tuating avalanches is N^(s)dX = NPx(s)(s + l)^r- Therefore, the conditional probability 
distribution of A backward avalanches obeys: 



and 1 @. 



is given by (s + Increasing the parameter to A + dX, the number N of avalanches 



d\n(s) x 
dX 



(£+1) 
1-A ' 



(5.1) 




c(s + l)P x (s) 
(1-A) 



c 



d(X 



2 Px(s)) 



(5.2) 



A(l -2A) 



dX 
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The proportionality constant c is determined from the normalization condition Y^^Lo P\( s ) = 
1, so that c = 

The distribution of all backward avalanches measures avalanches which begin at every 
time step s', and end at the first moment backward in time, s' — s, when \ m in{ s ' — s) > 
0'). We find 



P 



all i 



r^= x) ^ s)dX = 2 S v\m + i) + sprnj- (5 - 3) 



At large times this distribution is a power law Pfis) ~ s~ T Z where rf = 3/2. Here 
we have explicitly demonstrated that time reversal symmetry is broken. The forward and 
backward avalanches have different probability distributions and different scaling limits at 
large times. Eq. ( |5.3|) is in perfect agreement with numerical simulations. 

Unless the backward avalanche has size s = 1, two extremal values that are chosen at 
subsequent times must have both been present when the first of them was chosen. These two 
barrier values will have some ultrametric distance between them. This ultrametric distance 
must be less than or equal to the backward avalanche size, because the extremal value 
that terminates a backward avalanche must be an ancester to all of the extremal values in 
that backward avalanche. The probability to have an ultrametric distance larger than s is 
bounded above by P^ ll (s). We have numerically measured the ultrametric distance between 



subsequent activity and find a power law. In fact, as shown in Fig. [12] it has the same 
leading asymptotic behavior as P^ u within numerical accuracy. 

We conclude by speculating about connections with other phenomena related to glassy 
dynamics (see also Ref. |2T|. The Directed Polymer in a Random Media (DPRM) ||29|| , which 



could be a paradigm for glassy systems, exhibits an ultrametric structure in the optimal paths 
as well as a non-Gaussian tail for the probability distribution of these paths. These paths are 
somewhat analogous to the activity pattern in our model. However, unlike the DPRM, our 



model is inherently dynamical. Tang and Bak |30| found stretched exponential relaxation 



for the current in a sandpile model of SOC which could indicate that glassy dynamics takes 



place near a critical point. Recently Stein and Newman have put forward a picture of 
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dynamics on a high dimensional rugged fitness landscape based on an invasion percolation 
(SOC) picture. 

Finally, we note our model may fit into the picture of hierarchically constrained dynamics 
put forward by Palmer, Stein, Abrahams, and Anderson [|32] . We have an equation of motion 
for the dynamics which takes place in terms of avalanches spanning all time scales. These 



one 



avalanches are our hierarchically constrained degrees of freedom. Looking at Fig. |10 
notices that every A avalanche is composed of sub-avalanches which are fully contained 
within it. Each A avalanche cannot terminate until its sub-avalanches finish, so that the 
faster degrees of freedom successively constrain slower ones and form a hierarchy. 
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APPENDIX A: SOLUTION OF EQ. Q 

It is straightforward to solve Eq. (|3.3|) using a generating function. Defining p\(x) = 
J2^Lo xS P\( s ) an d applying the initial conditions [see Eq. (|3~2"1)1 P\(0) = and P\(l) = 
(1 - A) 2 , we find 

Px(x) - (1 - A) 2 = 2A(1 - \)xp x (x) + \ 2 xp x {x)\ (Al) 
an ordinary quadratic equation for the generating function. Its only acceptable solution is 

2 



(l - yjl -4A(1 - \)x 



Px(x) = ± '-. (A2) 

The solution for Px{s) in Eq. (|2.2j ) is simply given by the coefficients in x of the Taylor series 
of px{x). 
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The generating function p\(x) has a square- root singularity which determines the asymp- 
totic behavior of its Taylor coefficients, i. e. P\{s). We point out that this asymptotic 
behavior is a robust feature with respect of changes in the way the model is updated at 
each time step. For instance, if we had not chosen to set the barrier with the current min- 
imum to unity but to replace it also with a new random number, we would have obtained 
a cubic equation for p\(x); the solution of which would still be dominated by a square-root 
singularity. Furthermore, an update including more than nearest-neighbor sites would lead 
to even higher-order algebraic equations for p\(x), which are still dominated by the same 
square-root singularity. It would be interesting to consider changes in the updating rules 
that in fact would replace the leading square-root singularity, and the physics that such new 
rules indicate. 

APPENDIX B: CALCULATION OF EQ. O 

Let r be a nonnegative integer. Let N(i, r) be the probability that a A c avalanche which 
starts at site % will always be completely contained inside a box of radius r centered at the 
origin, % — 0. By definition, N(i,0) = because in the initial state the avalanche already 
contains one barrier which certainly will never be contained inside a box of vanishing radius. 

The properties of an avalanche starting at i can be related to the properties of avalanches 
starting at % — 1 and % + 1 by considering all possible states of the avalanche after the first 
step. No new active barriers are created after the first step with probability (1 — A c ) 2 = 1/4 
and the avalanche terminates without ever spreading beyond the site i. A single new active 
barrier is created after on time step either to the left or the right of i, each with a probability 
of A c (l — A c ) = 1/4, and N(i, r) is related to N(i — 1, r) or N(i + 1, r), respectively. Finally, 
with a probability of A 2 = 1/4 two new active barriers are created at i — 1 and at i + and 
the avalanche starting at i will never leave the box if neither avalanche ensuing from i — 1 
and i + 1 will ever leave the box. Thus, assuming throughout that r > 1, it is for |z| < r — 1 

N(i, r) = - + - [N(i - 1, r) + N(i + 1, r)} + ^N(i - l,r)N(i + 1, r). (Bl) 
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Clearly, N(i,r) = for all |z| > r, leading to the boundary condition for % = r — 1, 

N(r-l,r) = - + -N(r-2,r). (B2) 

Since N(i,r) is symmetric in its first argument, we will only consider nonnegative values of 
i and obtain another boundary condition at i — from Eq. (|B 1\) : 

N(0, r) = \ + \n{1, r) + ijV(l, r) 2 . (B3) 



We can simplify Eq. ( |B1|) by substituting 

jV(i,r) = l-/(i,r) (B4) 

to obtain for < i < r — 1 

A 2 /(z,r) = ~/(z-l,r)/(z + l,r), (B5) 
where Aj is a difference operator. Eq. (|B3Q leads to a boundary condition at i — 0, 

A,/(i = l,r) = i/(l,r) a , (B6) 
while Eq. (|B2|) gives another boundary condition for i — r — 1, 



Ai/(i = r - 1, r) = -^/(r - 2, r) + i. (B7) 

To make further progress we assume r>l, which allows to consider the continuum limit 
of Eq. (p5|). We set z = i/r such that z is a continuous variable in the unit interval for any 
r in this limit. We also set y(z) = f(i, r), where y is a function of z that depends on r as a 
parameter. Thus, we can rewrite Eq. ([Bo]) as 



- 2 y"{z) = \y{zf (B8) 



to leading order in the limit r — > oo. For the boundary conditions at z = and z = 1 we 
find from Eqs. (|B6^[B7D 



V(0) = ^(0) 2 , V(l) = -^(l) + i. (B9) 
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We note that we can obtain an equation for any value of A: 



2)y{z) + Xy(zf 



(BIO) 



which reduces to Eq. ( |B8| ) for A = A c . It is easy to show that the linear term dominates 
on the right-hand side for A < A c , leading to avalanche distributions with an exponential 
cut-off for large r. 

We can integrate Eq. (|B8|) using standard techniques for autonomous ordinary differential 



equations |25|| . We set u(y) = y' r (z), use the chain rule to get y"(z) = u'(y)u(y), and integrate 
once to find 



-y\z) = ± x -y(zf + C. 



(Bll) 



Since y(z) is a rising function of z, we choose the positive root. The integration constant C 
can be rewritten using the boundary condition at z = in Eq. (|B9|) as C = — |y(0) 3 , where 
we neglected terms of higher order in y(0) because y(0) is expected to be small for r — > oo. 
Integrating one more time we obtain 

d( 



y( z ) 
«/(") 



zr\ 



'y(o) 



(B12) 



"V 3 

We find at z — 1, using Eq. flB"9|), that y(l) = 2/3, because y'{l)/r can be shown to be 
negligible. Thus, y(l)/y(0) — > oo, and we obtain Eq. 



/(0,r)~y(0) 



3 / r°° d( N 






2 1 17.69 






[ r(f) . 


^2 jyi2 



(B13) 



Note that this result verifies our assumption that y(0) is small. Using dominate balance 
techniques |25] we can show that this solution is in fact the only consistent solution. 



APPENDIX C: LEADING ASYMPTOTIC BEHAVIOR OF THE 
SPATIO-TEMPORAL CORRELATIONS WITH RESPECT TO A PARTICULAR 

SITE 

The nonlinear integro-difference equation. ( |3.7| ) can be solved exactly in the continuum 
limit. Taking the continuum limit is justified because we are ultimately interested in the 
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behavior of F(r, s) for sufficiently large values of r and s. In general, to obtain the asymptotic 
behavior of a difference equation from the corresponding differential equation can be tricky 
even in the linear case P5] . But comparison with our calculation in Appendix D, where the 
continuum limit poses no problem, independently confirms our approach here. 

With f(r,x) = Yl^Lo x s F\(r, s ) an d p(x) = J2^Lo xS ^( s ) as generating functions, we 
obtain from Eq. ( |3.7|) , using Eq. (|A2"D , 



A 2 J(r, x) = A{x)f{r, x) + B{x)f{r - I, x)f{r + l,x), 

f(0,x)=p(x), and /(oo,x)=0, (CI) 



defining 



!-2xA(l-A)-2xAVx) _ x\> 

y ) x\{l - X) + x\ 2 p{x) ' y ) xX{l - X) + xX 2 p{x)' y ' 



We can take the continuum limit of Eq. ( |Cl| ) and get to leading order for large r an 



ordinary second-order nonlinear differential equation for / as a function of r: 

f(r)» = Af(r) + Bf{rf, /(0) = p, /(oo) = 0, (C3) 

where we have suppressed dependence on the parameter x. Using again the techniques for 
autonomous equations (see Appendix B) and the fact that /(oo) = 0, Eq. (|C3[) can be solved 
exactly to give 



f(r) =p 



cosh ( -\/~Ar \ + </l H — sinh ( — y/Ar) 

\2 J V V2 / 



(C4) 



Thus, we obtain from Eq. (|C4|) a closed-form expression for the envelope function of the 
spatio-temporal distribution of avalanches: 

F A (r )S )~-/^s— Vfas), (C5) 

where the contour encircles a small neighborhood of the origin in the complex x-plane in 
the positive direction. 
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From now on, we only consider the critical case, A = A c = 1/2. The integrand in Eq. (|C5[) 
can be expanded for large s in the neighborhood of x = 1 by substituting x = 1 — u/s. Then, 
the integration for u follows a contour that crosses the positive real axis from above near 
the origin in the complex-w plane. With 



fu u . fu „it „ 1 1 fu 



we find 

F(r, s) - 

where an analysis of Eq. (|C4| ) yields 



c 2-nis ^ ( s J 



(C6) 



(C7) 



12 



U , 



f(r, 1 - -) ~ 

5 



u \ 2 / r 



it \ 2 r 



h 

756 ~ 
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U \ 2 



exp 



(f ) 



1 < r 4 < s), 



(C8) 



In the first case of Eq. (|U8|) we have neglected terms with integer powers in u which would 
vanishes in the following integration. 

For r 4 ^> s, we evaluate the integral for F(r, s) by steepest-descent analysis similar 
to Sec. Ill, but taking account also of the non- exponential factor in the integrand. For 
1 C r 4 C s, we use Hankel's contour integral representation of the T-function |33] 

1 f du 



-v) 



c lui 



-e u u v 



(C9) 



to find 



F(r, s) 



J_ s -| fiii i J_rl _l 

v^F 6 \ L + r 2 + 84 a + 



f^^)"^ -1(4 



(K r 4 < s), 
fr 4 > s > 1), 



(CIO) 



where the second case is our main result for the non-Gaussian tail given in Eq. ( |3.9| ). 



APPENDIX D: LEADING ASYMPTOTIC BEHAVIOR OF THE 
SPATIO-TEMPORAL CORRELATIONS WITH RESPECT TO A BOX 



In this Appendix we calculate the distribution for the size of a box in space and time 
such that a A- avalanche will be fully contained in it. This notion, which easily generalizes 
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to higher dimensions, extends on the calculation in Sec. Ill and Appendix C, where we only 
considered avalanches with respect to a single site. While the calculation here (which is 
similar to Appendix B) is somewhat more extensive, the results are virtually identical and 
justify our simplified approach in Sec. III. 

Let P(r, s, i) be the probability for a A-avalanche, which starts at time s = with a 
single active barrier on site i, to have no active barriers for the first time at time s, and to 
not have left a box i G (— r, r). We merely consider avalanches in the critical state and set 
A = A c = 1/2. By definition, P(r = 0, s, i) = for all s > 0. Furthermore, P(r, s = 0, i) = 
for all r, because by definition no avalanche ends at time s = 0. Ultimately, we want to find 
the distribution F(r, s, 0) of avalanches of duration s that start at the origin i = and are 
completely contained in a box of radius r. A generic avalanche is plotted in Fig. |J The 
smallest box it is fully contained in is of size r = 18 in this case. In correspondence with 
Eq. (fUD it is 

F(r, s, i) = P{r = oo, s, i) — P(r, s, i), (Dl) 

where P{r = oo, s, 0) = p(x) is given in Eq. (|A2|) . 

As before, the properties of an avalanche that originates at s = can be deduced from the 
properties of avalanches that ensue after the first update. The original avalanche can either 
terminate after the first update when the update does not produce any new active barriers 
with probability (1 — A) 2 = 1/4, or it can generate new avalanches by creating new active 
barriers. If the first update creates exactly one new barrier with likelihood A(l — A) = 1/4 
either to the left or to the right of site i, the properties of an original avalanche of duration s 
is related to the properties of an avalanche of duration s — 1 with regard to a site of distance 
i— 1 or i+ 1, respectively. If the first update creates two new active barriers with probability 
A 2 = 1/4 to the left and the right of site i, two new avalanches ensue. Then, the properties 
of the original avalanche of duration s is related to the properties of all combinations of 
two avalanches of combined duration s — 1. For any such combination, the probability to 
not leave the box when starting at site i is given simply by the product of the probabilities 
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for the two ensuing avalanches to not leave the box after starting at site % — 1 or % + 1, 
respectively. We thus obtain for r > 1 and \i\ < r that P(r, s — — 1/4, and for all s > 1 
that 

P(r, s + 1, z) = - [P(r, s, i - 1) + P(r, s, z + 1)] 

+~ £ P(r, a', i - l)P(r, s - s', z + 1). (D2) 

4 s'=0 

Since P(r, s,i) is symmetric in z, we restrict ourselves to nonnegative values of z, leading to 
a boundary condition at i = 0: 

P(r, s + 1, 0) = ~P(r, s,l) + \j2 P(r, s', l)P(r, s - s', 1) (s > 1, r > 1). (D3) 

s'=0 

Since P(r, s, z > r) = 0, we obtain a second boundary condition at i — r — 1: 

P(r,s- l,r- 1) = ^P(r,s,r-2) (s>l,r>2). (D4) 
Using Eq. (|D1|), and defining /(r, x, z) = Y^Lo x s F(r, s, z), we obtain for 1 < i < r — 2 
A 2 /(r, x, z) = A(x)f(r, x, z) + B(x)f(r, x, i - l)/(r, x, z + 1) (D5) 
where A(x) and P(x) are the same as in Eq. (|C2] ), supplemented by the boundary conditions 



Aif(r,x,i = 1) 



3j 3j / > 

1 p(x 

4 4^ ; 



/(r,x, 1) + -f(r,x, if 



Ai/(r, x, i = r - 1) = (| - l) /(r, x, r - 2) - p{x) (l - | ) - |. (D6) 

As before in Appendix C, we expand the equations in the limit x — > 1 to analyze the 
avalanche distribution for large times s. Then Eq. ( |D5|) simplifies for 1 < z < r — 2 to 

1 



A 2 /(r, x, z) - 2 v / I^/(r, x, i) + -f(r, x, i - l)/(r, x, z + 1) (D7) 

with the boundary conditions 

Ai/(r, x, z = 1) ~ vT^/(r, x, 1) + i/(r, x, l) 2 , 
Ai/(r, x, z = r - 1) ~ -~f{r, x, r - 2) + ~. (D8) 
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For sufficiently large r, we can take the continuum limit of these equations where f(r, x, i) — > 
y{z)jr 2 with z = i/r as a continuous variable in the unit interval. Eq. ( |L)7| ) then approaches 



y"{z) = 2r 2 VT^y(z) + ^y(z) 2 , y'(0) = rVT=^/(0), y{l) = jjr 2 . 



(D9) 



We can obtain a first integral of Eq. ( |D9j) using again the technique for autonomous equa- 
tions: 



y'(z) = | |j,(z)3 - y( )3] + 2r=(l - x)* - y(0) 2 ], 



(D10) 



assuming that j/'(0) — > for r — ► oo sufficiently fast. Since y(z) is a rising function of z, we 
have to chose the positive root. Integrating again, and using y(l), we obtain 

6r 2 \/l — x 



.%( 0) 



'v(o) 



/i ^e-i + ^^-i) v 3 

For y(0) <C r 2 vT^~x, i. e. a > 1, we get 

1 



z with a 



y(o) 



(Dll) 



'y(o) 



~ — p= In 



3 

2* 



(D12) 



which yields 



f(r,x,0) 



2/(0) 



3^ — x exp (-V2r(l-x)*) (K(l 



(D13) 



and by steepest-descent analysis as before 

F(r, s, 0) ~ C(r, s) exp 



(D14) 



confirming the non-Gaussian tail found in Eq. ( ^.9| ). A similar consideration of the integral 
in Eq. (|D 1 1| ) would determine the behavior for r 2 ^> y(0) ^> r 2 \/l — x, i. e. a <C 1. 
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FIGURES 

FIG. 1. Punctuated equilibrium behavior for the evolution of a single species in the 
one-dimensional M = oo model. The vertical axis is the total number of returns of the activ- 
ity to some site as a function of time s. Note the presence of plateaus (periods of stasis) of all 
sizes. The distribution of plateau sizes scales as s~ 7 / 4 . 
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FIG. 2. Snapshot of the stationary state in a finite one dimensional system for the M = 1 
(Bak-Sneppen) model. Except for the avalanche which consists of small fitness values in a localized 
region, almost all the fitness values in the system are above a self-organized threshold A c . 
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FIG. 3. Ultrametric tree structure. At any given time, indicated by the vertical axis, all of the 
active sites below threshold have an ancestry which forms a tree. The ultrametric distance between 
any pair is the distance back in time to the first common ancestor. 
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FIG. 4. Numerical verification of Eq. (2.4). We sampled the probability for an avalanche to 
start at the origin and to be contained by the smallest possible box of radius r for r < 250 with 
increments of Ar = 10. Squares represent the measured values from a sample of 2 x 10 6 avalanches 
for this probability, corresponding to —d r f(0, r) ~ 35.48 . . . r~ 3 which is drawn as a solid line. 
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FIG. 5. Plot of a typical A c avalanche for M = oo, starting at the origin at time s = 0, 
and ending at time s = 976. At each time step, every site with at least one active barrier is 
marked (}, and the sequence of minimal sites is connected with a line showing jumps of various 
sizes. Apparent also are the punctuated equilibria for each site which extend over many sizes. 
The propagator G(r, s) of the activity is the probability to have a minimum at site r at time s, 
while F(r, s) is the probability for an avalanche to end at s (here: = 976) and to have reached a 
particular site r during its duration. 
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FIG. 6. Counting rates for the activity plotted as a function of time s. Each curve is labeled 
in descending order by i = [log 2 x\ where x = r 4 /s is the scaling variable [see Eq. ( |3.17|) 1. Data 
with i > 12 has been disregarded due to insufficient statistical accuracy. 
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FIG. 7. The normalized distribution of activity G(r, s) derived from Fig. ||. Each value of 
G(r, s) is divided by its corresponding value of r and cut off below s = 10 4 and above s = 5x 10 . 
All data "collapses" onto plateaus with a mean value for each plateau given by x on the right, 
again labeled by i = [log 2 x] . The standard deviation for the data in each plateau is indicated by 
error bars. 
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FIG. 8. Plot of the extrapolants from Eq. ( 3.18Q for each i = |4og 2 x\ as function of 1/lnx. 
Asymptotically for 1/lnx — > 0, this sequence approaches the value of a. From the extrapolation 
of the sequence to x — > oo (indicated by the lines that extend the sequence and its errors to the 
ordinate) we estimate a = 0.35 ± 0.03, in good agreement with a = 1/3. 
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FIG. 9. Plot of the mean width R(s) of the compact region of covered sites for surviving 
avalanches at time s for s < 10 5 . The result of a simulation involving 10 s avalanches is given by 
0, and the solid line is given by 6.10 (^s^ — 1.0^, obtained from an asymptotic fit of the data for 
large s. 
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FIG. 10. Plot of a part of the sequence of minimal random numbers A m in(s) chosen for an 
update at time s in a A c avalanche for M = oo. The durations of a hierarchy of A avalanches is 
indicated by forward arrows, where A = A m i n (s). Longer avalanches with larger values of A contain 
many shorter avalanches which have to finish before the longer avalanche can terminate. Note that 
an update punctuates any A avalanche with A < A m i n (s). 
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FIG. 11. Same sequence A m i n (s) as in Fig. |T^, where the durations of backward A avalanches is 
indicated by backward arrows. A similar hierarchical structure of sub-avalanches emerges, although 
with a different distribution than for forward A avalanches in reflection of the irreversibility of the 
update process. 
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FIG. 12. Log-log plot of the distribution for the duration of backward avalanches (+) and the 
distribution for the ultrametric distances between subsequent minimal sites (x). Both functions 
seem to coincide asymptotically for large s. 




